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We estimate, using a large-scale Monte Carlo simulation, the critical expo- 
nents of the antiferromagnetic Heisenberg model on a stacked triangular lat- 
tice. We obtain the following estimates: 7/1/ = 2.011 ± .014, v = .585 ± .009. 
These results contradict a perturbative 2 + e Renormalization Group calcu- 
lation that points to Wilson-Fisher 0(4) behaviour. While these results may 
be coherent with 4 — e results from Landau- Ginzburg analysis, they show the 
existence of an unexpectedly rich structure of the Renormalization Group 
flow as a function of the dimensionality and the number of components of 
the order parameter. 
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1 INTRODUCTION 



There is at present a satisfactory understanding of the critical behaviour of 
physical systems where the rotation symmetry group O(N) is broken down 
to O(N-l) at low temperatures. Several theoretical tools are available to 
estimate the critical exponents and there is good agreement between these 
estimates. The Wilson-Fisher fixed point which describes the critical physics 
can be smoothly followed between two and four dimensions: for iV > 3 the 
2+e and 4— e renormalization group expansion merge in a continuous manner. 

The situation is much more complicated when the rotation symmetry is 
fully broken in the low-temperature phase. A prominent example is found 
in the so-called helimagnetic systems where Heisenberg spins are in a spiral 
arrangement below the critical temperature. It is an interesting question, 
both theoretically and experimentally, to know the corresponding univer- 
sality class. A widely studied prototypical model is the antiferromagnetic 
Heisenberg model on a stacked triangular lattice, which is simple and dis- 
play commensurate helimagnetic order below a transition point T c . There 
is little consensus in the literature regarding critical phenomena associated 
with this model 0. 

This topic has been investigated by use of a D = 4 — e renormalization 
group calculation^. The corresponding Ginzburg-Landau theory for a N- 
component vector model involves two N-component bosonic fields. It is found 
that for N large enough the transition is second order and not governed 
by the Heisenberg 0(2N) Wilson-Fisher fixed point but by a different fixed 
point which is also non-trivial for D < 4. For smaller N, this new fixed 
point disappears and there is no stable fixed point which is an indication 
for a fluctuation-induced first-order transition. The dividing universal line 
between second-order and first-order behaviour is found to be N C (D) = 21.8— 
23.4e + 0(e 2 ). The rapid variation of N c leaves us rather uncertain about the 
fate of the case D = 3. Clearly more information is needed about the N C (D) 
line in the (number of components-dimension)-plane. 

AD = 2 + e renormalization group study has been performed for a 
system of Heisenberg spins by use of a non-linear sigma model defined on a 
homogeneous non-symmetric coset space 0(3) x 0(2) jO (2). It was found that 
near two dimensions the system undergoes a second order transition which is 
governed by the 0(4) usual Wilson-Fisher fixed point. In fact the symmetry 
0(3) x 0(2)/0(2) is dynamically enlarged at the critical point to 0(3) x 
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0(3)/ 0(3) and 0(3) x 0(3) is 0(4). This mismatch between the expansions 
near four dimensions and near two dimensions is quite unusual and does not 
happen in the well-studied 0(N) — > 0(N — 1) critical phenomena. As a 
consequence, the D = 3 case remains elusive and a direct study in three 
dimensions is called for. 

Some preliminary Monte Carlo (MC) simulations have shown evidence 
[|], [5], |6|, 0] for a continuous transition in the case of the Heisenberg anti- 
ferromagnet on a stacked triangular lattice. The exponents found are not 
compatible with those of the 0(4) vector model in three dimensions. We did 
a large-scale simulation with much better statistics than previous attempts. 
As a consequence we are able to pin down the transition temperature in a 
very precise manner and to obtain reliable estimates of critical exponents. 

We thus focus on the classical spin model defined by the classical Heisen- 
berg Hamiltonian: 

H = J ij S » " S i- (!) 

<ij> 

The exchange interaction is nonzero ( Jy = 1 in what follows) between 
nearest-neighbors of a stacked triangular lattice and the spins are three- 
component unit vectors. The classical ground state is found by minimizing 
the Fourier transform J(Q) of the couplings Jy. The spins adopt a planar 
arrangement on a three-sublattice structure with relative angles 120 degrees. 

2 The simulation 

We made two successive sets of Monte Carlo simulations of this model, using 
slightly asymmetric lattices (with periodic boundary conditions) of shape 
L 2 L Z , where L is the linear size inside the planes, and L z = 2/3L the stacking 
size. 

The first set of simulations was run on a CM-2 8K massively parallel 
computer. We made runs on lattices of size 12 2 8 (with a total of 600 000 
Monte Carlo sweeps of the lattice), 24 2 16 (with a total of 1 500 000 sweeps) 
and 48 2 32 (with a total of 1 570 000 sweeps). In the first case we used 
the Heat-Bath algorithm, in the later two cases we used an Hybrid Over- 
relaxation algorithmic where each Heat-Bath sweep is followed by an energy 
conserving sweep. The next set of simulations was run on CRAY vector 
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computers on lattices of size 18 2 12, 24 2 16, 30 2 20 , 36 2 24 and 48 2 32 with 4 
to 8 10 6 Hybrid Overrelaxation sweeps. Both simulations concentrate in the 
immediate vicinity of the transition. 

Let us note by S a the total spin per site for sublattice a (a £ [1, 2, 3]), we 
measured the magnetization: 

M L ([3) = W<\Sa\>, (2) 

a 

the susceptibility: 

xM = ^ (£<s 2 >-<|S a |> 2 ), (3 ) 

and the fourth-order cumulant: 

Our strategy to extract the critical exponents was first to estimate the 
value of the critical temperature /3 C , as the point where Bl({3) is L indepen- 
dent and then to estimate the exponents from the following set of equations: 



XlWc 
M L ((3 C 
dB L {(3) 



8(3 

|>nM l(/3 ) 



P=Pc 



L i/u 



(5) 
(6) 

(7) 
(8) 



The extrapolation, from the j3 value used for the simulation, to /3 C is done 
using the Ferrenberg-Swendsen technique |TIJ. This technique is invaluable 
to extrapolate in a neighborhood of size ~ l/L 1 ^ (We understand that, 
when used blindly outside such a tight range, it may give wrong results). 
An alternative strategy would be to evaluate exponents from the finite size 
scaling of the maxima of Xl((3), M l (/3), -j^B L (f3) and -^\nM L (f3). It would 
however require to perform simulations at four (3 values, close to the four 
(L dependent) points where this quantities reach their maxima, otherwise 
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accuracy would be lost. The statistical analysis is done with respect to 20 
bins, using first-order bias-corrected jackknife (see e.g. Ref. ]TI| and references 



therein). The first 20 % of each run is discarded for thermalization. 

From simulations of lattices of increasing sizes L\ < L2 < ■ ■ ■ converging 
estimators of (3 C are obtained by solving the equations: 

13^0) = B Li 0). (9) 





/3 C 


j/v 


D - 2[3/v 


V 


V 


24-12 
48-24 


1.04401 (38) 
1.04462 (21) 


2.009 (14) (04) 
2.033 (18) (14) 


2.030 (07) (07) 
1.998 (09) (25) 


.651 (22) (03) 
.590 (20) (04) 


.604 (10) (01) 
.576 (10) (04) 


24-18 
30-24 
36-30 
42-36 
48-42 


1.04452 (40) 
1.04388 (36) 
1.04408 (37) 
1.04485 (25) 
1.04427 (30) 


2.010 (26) (06) 
1.994 (36) (11) 
1.988 (55) (16) 
2.083 (68) (14) 
1.986 (61) (18) 


2.011 (08) (10) 
2.021 (13) (16) 
2.020 (18) (21) 
1.984 (20) (31) 
2.004 (22) (39) 


.605 (30) (03) 
.578 (35) (03) 
.605 (55) (12) 
.528 (48) (06) 
.639 (69) (05) 


.608 (17) (02) 
.583 (19) (03) 
.607 (30) (04) 
.548 (28) (06) 
.577 (27) (03) 


48-18 


1.04432 (06) 


2.011 (07) (12) 


2.010 (03) (21) 


.590 (09) (06) 


.588 (04) (03) 


all 


N/A 


2.011 (07) (12) 


2.011 (02) (19) 


.585 (07) (06) 


.585 (04) (04) 



Table 1: Critical temperature and critical exponents estimates, using equa- 
tions 9, 5, 6, 7 and 8 respectively. First two lines are CM-2 data, next lines 
the CRAY data. The first column gives the linear sizes of the lattices used. 
The two numbers inside parenthesis are the estimated statistical errors, direct 
and induced, on the last digits of the number on their left. 

The results can be found in the second column of Tab.l. Within our 
statistical accuracy all estimates are compatible, there is neither any clear 
lattice size dependence, nor any discrepancy between the two sets of runs. 
Our final estimate is (the quoted error is deliberately on the safe side) 

p c = 1.0443 ± .0002. (10) 

We use this value to compute Xl(Pc): M l ((3 c ), -^B l (/3) and -^\n.M L {(3) 
in order to estimate 7/V, j3/v and v . Results can be found in Tab.l. We 
quote separately the estimated "direct" statistical errors computed from the 
dispersion of the results from the 20 bins, and the errors induced by the 
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uncertainty on the determination of /3 C , computed as the (absolute value of 
the) difference between the values obtained using our best estimate of f3 c and 
the values using the one standard deviation estimate from Eq.|lO|. The last 
line in Tab.l gives the results of linear fits of the CRAY data for \ii(xl(Pc)) , 

ln(M L (/3 c )), H^B L ((3) 



) and ln(jUnM L (/3) 



) respectively. The re- 

suits of the fits are stable against omitting the smallest lattice data. The x 2 
is equal to 1.5, 3.7, 1.9 and 4.5 respectively (we expect x 2 ~ 4). Such values 
mean that linear fits give good representation of the data for L > 18 (with 
the current accuracy). 

Our final numbers are 7/1/ = 2.011 ± .014, D - 2f3/v = 2.011 ± .019, 
v = .585 ± .009. We observe that hyperscaling is verified within errors, and 
1] is very small (jr/l ~ .01). The value for v is clearly not compatible with the 
value for the 0(4) fixed point Hl3l [y ~ .75). 

Published Monte Carlo results for v are .53±.02||, .53±.03@, .55±.03§ 
and .59 ± .020. Our results are in good agreement with the most recent re- 
sults (and with the result .57±.02 from a model of commensurate Heisenberg 
helimagnet||]). The methods of analysis are quite different. Ref.0 uses data 
taken in a wide region around the transition point and adjust the values of 
the critical exponents using the old fashioned "data collapsing" method. This 
method has the disadvantage to give much weight to points with large value 
of (P — (3 c )L l l u '. We extract the exponents directly from data reweighted 
to (3 C . We have also a much higher statistics: we use 8 x 10 6 heat-bath + 
energy-conserving sweeps whereas uses 6-20 times 20000 sweeps of a less 
efficient algorithm. Our estimated statistical error on the determination of 
(3 C is one order of magnitude smaller than the one in ref. W\. 



3 Conclusion 

Let us first summarize the knowledge obtained from perturbative RG studies 
as a function of the dimensionality and the number of components of the order 
parameter. 

i) There is a universal line N C (D) separating a first-order region from a 
second-order region near D = 4 — e. Its slope is known from RG studiesfQ] 
as well as the critical value iV c for D = 4. 

ii) Large-N studies[fr2[ indicate that the stable fixed point found above N C (D) 
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persists smoothly in the region N = oo and 2 < D < 4. 
iii) Smoothness along the D = 2 vertical axis can be shown by studying 
the nonlinear sigma model suited to the N-vector model: it is built on the 
homogeneous space 0(N) x 0(2) /0(N — 2) x 0(2)0. One finds a single 
stable fixed point for all values of iV > 3. In the large-N limit the exponents 
from this sigma model are the same as those of the linear model. 

If we believe in the perturbative RG results, then necessarily the universal 
line N C (D) can only intersect the horizontal axis N = 3 between D=2 and 
D=4. The simplest hypothesis is then that the plane (N, D) is divided in two 
regions by the line N C (D). This line may intersect the N = 3 axis at a critical 
dimension D c . Our Monte-Carlo results thus imply that D c is between three 
and four dimensions since we observe a continuous transition at D=3. A 
possibility would be that the critical behaviour we observe is described by 
the fixed point found in the neighborhood of D=4. 

However one has to note that the sigma model approach shows that this 
fixed point is 0(4) for N=3. This is clearly ruled out by our data. The 
simplest scenario suggested by i-iii above is wrong. It may be that the whole 
sigma model approach breaks down since perturbative treatment of a sigma 
model neglects global aspects. Previous studies indeed have performed, as 
usual, only perturbation theory for spin-wave excitations^]. One may invoke 
the peculiar vortices present in the model Jl6j since Ui(S0(3)) = Z2 that 
form loops in D=3 as responsible of the breakdown of the sigma model. 
However vortices are also excluded from the perturbative calculations near 
D=4 and thus it would be difficult in this case to believe in what is found 
perturbatively. It may be also that the fixed point found in 2 + e expansion 
is destabilized by some operators containing high powers of gradients [IS|. 

We have thus obtained the critical behaviour of a magnet with a canted 
ground-state. This calls for a deeper understanding of the Renormalization 
Group behaviour of canted systems since there is no obvious relationship be- 
tween the different perturbative approaches, contrary to the case of collinear 
vector magnets. We note for the future that it would be interesting to obtain 
the critical behaviour of the XY canted systems, much more relevant to the 
experimental situation. 



7 



Acknowledgements 

We would like to thank Guy Decaudain, David Lloyd Owen, Gyan Bhanot 
and Roch Bourbonnais of TMC for their support. We also thank H. T. Diep, 
A. P. Young and D. P. Belanger for various informations about this subject. 

References 

[1] For a review see: M. L. Plumer and A. Caille, J. Appl. Phys. 70, 5961 (1991). 

[2] T. Garel and P. Pfeuty, J. Phys. C9, L245 (1976); D. Bailin, A. Love and M. 
A. Moore, J. Phys. CIO, 1159 (1977). 

[3] P. Azaria, B. Delamotte and Th. Jolicoeur, Phys. Rev. Lett. 64, 3175 (1990). 

[4] H. Kawamura, J. Phys. Soc. Jpn 54, 3220 (1985). 

[5] H. Kawamura, J. Phys. Soc. Jpn 56, 474 (1986). 

[6] H. Kawamura, J. Phys. Soc. Jpn 58, 584 (1989). 

[7] H. Kawamura, J. Phys. Soc. Jpn 61, 1299 (1992). 

[8] H. T. Diep, Phys. Rev. B 39, 397 (1989). 

[9] A.D. Kennedy, Proceedings of the International Symposium on Lattice Field 
Theory LATTICE 92, September 15-19 1992 Amsterdam. 

[10] A.M. Ferrenberg and R.H. Swendsen, Phys. Rev. Lett. 61 2635 (1988); 

See also A.M. Ferrenberg and R.H. Swendsen, Phys. Rev. Lett. 63 1195 
(1989). 

[11] A. Billoire, R. Lacaze, and A. Morel, Nucl. Phys. B 370 773 (1992). 
[12] H. Kawamura, Phys. Rev. B38, 4916 (1988). 

[13] H. E. Stanley in "Phase transitions and critical Phenomena", Vol. Ill, C. 
Domb and M. S. Green eds., Academic Press, New York, 1974; J.C. Le Guil- 
lou and J. Zinn- Justin have also obtained unpublished estimates of the 0(4) 
exponents by field-theoretic methods that are very close to the series esti- 
mates. 



8 



[14] P. Azaria, B. Delamotte, F. Delduc and Th. Jolicoeur, Saclay preprint 
SPhT/93-044, to appear in Nuclear Physics B[FS]. 

[15] G. E. Castilla and S. Chakravarty, Phys. Rev. Lett. 71, 384 (1993). 

[16] H. Kawamura and S. Miyashita, J. Phys. Soc. Jpn 53, 4138 (1984). 



9 




DIMENSION 



